Discreteness-Induced Slow Relaxation in Reversible Catalytic 

Reaction Networks 



Akinori Awazu^ and Kunihiko Kaneko^'^ 
^Department of Mathematical and Life Sciences, 
Hiroshima University, Kagami-yama 1-3-1, 
Higashi- Hiroshima 739-8526, Japan. 
"^Department of Basic Science, University of Tokyo, 
Komaba, Meguro-ku, Tokyo 153-8902, Japan. and 
^ERATO Complex Systems Biology, JST, 
Komaba, Meguro-ku, Tokyo 153-8902, Japan. 
(Dated: March 1, 2010) 

Abstract 

Slowing down of the relaxation of the fluctuations around equilibrium is investigated both by 
stochastic simulations and by analysis of Master equation of reversible reaction networks consisting 
of resources and the corresponding products that work as catalysts. As the number of molecules 
N is decreased, the relaxation time to equilibrium is prolonged due to the deficiency of catalysts, 
as demonstrated by the amplification compared to that by the continuum limit. This amplification 
ratio of the relaxation time is represented by a scaling function ash = N exp(— and it becomes 
prominent as N becomes less than a critical value /i ~ 1, where /3 is the inverse temperature and 
V is the energy gap between a product and a resource. 

PACS numbers: 
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I. INTRODUCTION 



The study of reaction processes in catalytic reaction networks is generally important to 
understand the dynamics and fluctuations in biochemical systems and their functionality. 
Obviously, understanding the generic features of equilibrium characteristics and relaxation 
to equilibrium is the first step toward gaining such an understanding. Indeed, such reaction 
systems often exhibit anomalous slow relaxation to equilibrium due to some kinetic con- 
straints such as diffusion- infiuenced (limited) reaction [l] and formations of transient Turing 
patterns[2|. In this paper, we consider a novel mechanism to realize such slow relaxation in 
catalytic reaction networks, where the discreteness in molecule number that may reach zero 
induces drastic slowing down. 

Most intra-cellular reactions progress with the aid of catalysts (proteins), whereas cata- 
lysts have to be synthesized as a result of such catalytic reactions. Indeed, reaction dynamics 
in catalytic networks have been extensively investigated. In most such studies, a limiting case 
with a strong non-equilibrium condition was assumed by adopting a unidirectional reaction 
process (i.e., by neglecting backward reactions). To understand the basic properties of bio- 
chemical reactions, however, it is important to study both equilibrium and non-equilibrium 
characteristics by including forward and backward reactions that satisfy the detailed balance 
condition. Such a study is not only important for statistical thermodynamics but it also 
provides some insight on the regulation of synthesis or degradation reactions for homeostasis 
in cells. 

Recently, we discovered a slow relaxation process to equilibrium, which generally appears 
in such catalytic reaction networks, and proposed "chemical-net glass" as a novel class 
of nonequilibrium phenomena. In this case, relaxation in the vicinity of equilibrium is 
exponential, whereas far from it, much slower logarithmic relaxation with some bottlenecks 

n 

appears due to kinetic constraints in catalytic relationships p|]. In this study, we adopted 
continuous rate equations and assumed that the molecule number is sufficiently large. 

In biochemical reaction processes, however, some chemical species can play an important 
role at extremely low concentrations of even only a few molecules per cell ^-|6| . In such sys- 
tems, fiuctuations and discreteness in the molecule number are important. Indeed, recent 
studies by using a stochastic simulation of catalytic reaction networks have demonstrated 
that the smallness in the molecule number induces a drastic change with regard to statisti- 
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cal and spatiotemporal behaviors of molecule abundances from those obtained by the rate 
equation, i.e., at the limit of large molecule numbers j?- 20 1. In these studies, the strong 
nonequilibrium condition is assumed by taking a unidirectional reaction. 

Now, it is important to study the relaxation process to equilibrium by considering the 
smallness in the molecule number. Does the discreteness in molecule number influence the 
equilibrium and relaxation behaviors? Is the relaxation process slowed down by the smallness 
in the molecule number? To address this question, we have carried out several simulations 
of the relaxation dynamics of random catalytic reaction networks by using stochastic sim- 



ulations. Numerical results from several networks |21l. |22| suggest that the relaxation time 
is prolonged drastically when the number of molecules is smaller. The increase from the 
continuum limit is expressed by the factor exp{P6E), where 6E is the additional energy 
required to pass through the bottleneck due to the discreteness in molecule number and (3 
is the inverse temperature. 

In this paper, we analyze such slowing down of a reaction process to equilibrium that is 
induced by the smallness in molecule numbers. Instead of taking complex reaction networks, 
we choose simple networks or network motives to estimate the relaxation time analytically. 
In fact, complex networks are often constructed by combining a variety of simple network 
motives with simple branch or loop structures. We focus on the relaxation dynamics of 
reversible catalytic reaction systems with such simple network motives as a first step toward 
understanding the general relaxation properties in complex catalytic reaction networks. 

In section II, we introduce two network motives, where the synthesis of a product from 
resource molecules (and its reverse reaction) is catalyzed by one of the other products. Here, 
we note that some specific network motives may exhibit incomplete equilibration when the 
molecule number decreases, and the average chemical concentration in the steady state 
deviates from the equilibrium concentration derived by the continuous rate equations. 

In section III, we show relaxation characteristics from the stochastic simulations. The 
relaxation of the fluctuation around the steady state slows down as the molecule number 
is decreased below a critical value. This increase is represented by a scaling function by 
using h = N exp{—pV), where is the molecule number and V, the energy gap between a 
product and a resource. In section IV, we present an analytic estimate for this relaxation 
suppression due to the smallness in molecule number by using a suitable approximation for 
Master equation. In section V, we present a summary and discuss the generality of our 
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(a) Cascade 




FIG. 1: Illustration of (a) cascade system and (b) loop system. Solid arrows indicate reaction 
paths (their width indicates the transition tendency) and dashed arrows indicate catalyzation. 

results. 

II. MODELS 

Here, we consider reversible catalytic reaction systems with two simple network struc- 
tures, cascade system and loop system, as shown in Fig. 1, which may function as network 
motives for complex reaction networks. These systems consist of 25" chemical species, which 
are Product Pi and Resource Ri with i — 1,2, S. Here, each product chemical can catalyze 
at most one of the other Resource-Product reactions, whereas each reaction is catalyzed at 
most by some product. (Instead, we can interpret that there exist S chemical species with 
excited and non-excited states, and chemicals in an excited state can catalyze an excitation 
reaction of one of the other molecules.) 

If all chemicals are catalyzed by one of them, we can renumber and Ri for i — 
1,2,...,S — 1 and write the reaction as 

Pi + Pi+l ^kJ 1' Pi + Pi+li 

where Ps + Pi ^^^^'p^ Rs + Pi: which leads to the loop system (b). When there exist a 
reaction that is not catalyzed, the cascade system in Fig. la) is obtained where Ps ^fe^g p^ 
Rs- (Neglecting cases in which some pair of resource and product is totally disconnected 
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from others, the loop and cascade systems are the only possibilities). 

The rates of forward (kp^^R^) and backward (kn^^p^) reactions are set so that they satisfy 
the detailed balance condition. We assume that the energy of the chemical Pj is larger than 
that of Ri, and we set kp.^R^ = 1 and kji- p^ = exp(— /3Vi), where Vi is the energy gap between 
Pi and Ri and /3 is the inverse temperature. We define pi and rj as the number of molecules 
of the chemical species Pi and Ri, respectively. We fix the total number of molecules as SN, 
and pi + Vi — N holds for each i. The state of the system is represented by a set of numbers 
(pi,P2,-,Ps)- 

In both the systems, it is noted that for A?^ — >■ cxd (i.e., the continuous limit), < pi >— >■ 

Pi'^ = ^j^^-pv^ and < Ti >^ rl'^ = ^_^j^^y. holds at the equilibrium distribution, which is 
reached at t — )■ oo. 

For finite N , however, there is a difference between the distribution of the cascade and the 
loop systems. In the cascade system, the average of the equilibrium chemical concentrations 
are identical to the continuum limit, and are given by < pi >— y^^=Wi is, they are 
independent of N and ^. This is because all the states (pi,p2, ■■■,^5) (0 < Pi < are 
connected by reactions and the above equilibrium distribution is only the stationary solution 
for Master equation. 

On the other hand, in the loop system, there is a deviation in the steady chemical 
concentration from the continuum limit, which becomes more prominent as N becomes 
smaller. This is because the state (pi,P2, ■■■■,Ps) — (0, 0, 0) cannot be reached from other 
states, whereas the state cannot move to any other states. Hence, the steady distribution 
from the initial conditions without {pi,P2, ■■■,Ps) = (0, 0, 0) deviates from the continuum 
limit. This deviation becomes prominent as N becomes smaller. For example, for N = 1 
and Vi = V, the distribution from the initial condition without {pi,P2, ■■■,Ps) — (0,0, ...,0) 
is given by < pi >— ^ '^(^i^l-fiv'^sl^ — ■ Note that < Pi > tends to 1/S with an increase in (3. 

III. SIMULATION RESULTS 

In this section, we present the results of stochastic simulations and show the dependence 

of the relaxation process on the number of molecules and the inverse temperature /3. For 
simplicity, we consider Vi to be uniform for all species (= V); however, this assumption can 
be relaxed. 
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FIG. 2: C{t) of cascade systems with (a) 5 = 2, (b) 5 = 3, and (c) 8 = 4, and loop systems with 
(d) = 2 and (e) S = 3 for several N with /3 = 3. (f) C(oo) as a function of h in loop systems for 
several /3 and S. Code indicates the auto-correlations given by Eq. (4) in (a)-(c) and Eq. (3) in 
(d) and (e). C* = exp(-e-^^t) in (d), and C* = exp(-^^t) in (b) and (e) with /3 = 3{V = l). 

Numerical simulations are carried out by iterating the following stochastic processes, (i) 
We randomly pick up a pair of molecules, say, molecule 1 and 2. (ii) Molecule 1 is transformed 
with its reaction rate (if it is P, it is transformed to R, and vice versa) if molecule 2 can 
catalyze the reaction of molecule 1. In the cascade case, there is a reaction that progresses 
without a catalyst, and in this case, if molecule 1 is the one that reacts without a catalyst, 
then it is transformed with the reaction rate independently of 2. Here, a unit time is defined 
as the time span in which the above processes for catalytic reactions are repeated SN times. 
In each unit time, each molecule is picked up on average to check if the transformation 
occurs. 

In the following, we focus on the behavior of the system after a sufficiently long time 
from the initial time where the numbers of each molecule Pi and rj are set randomly from 
[0, A^] under the constraint pi + Vi = N and {pi,P2, ■■■,Ps) 7^ (0, 0, 0). 

Figure 2(a)-2(e) show the auto-correlation functions of the deviation from the equilibrium 
concentration of the cascade system ((a)-(c)) and the loop system ((d) and (e)) for some S 
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FIG. 3: r as a function of in (a) cascade system and (b) loop system, and t'^{S) and t^{S) for 
/3 = 3 and = 2, 3, 4. p as a function of h in (c) cascade system (d) loop system with S = 2,3,4 
for several /3 

and N with (3 = 3, defined by C{t) = c(t)/c(0) and c(t) =< Ei[fe(t) -^^^(0) - + 
(rj(t) — rj^'^)(rj(0) —r^'')] >■ As already discussed, C(oo) — > in the cascade system whereas 
C(oo) > for small A^. The value C(oo) starts to deviate when h = Ne~^^ becomes less 
than 1. Hence, we have plotted C(oo) of the loop system as a function of h in Fig.2(f) for 
/3 = 1 and 3. As shown, C(oo) > holds for /i < 1 independently of /3. On the other 
hand, in both systems, the relaxation to the final state with C{t) = const, for small is 
drastically slowed down as compared to that for large when 5* > 2, whereas the relaxation 
for small A^ is faster when S = 2. 

To observe the dependence of the relaxation time on A^, we measured the integrated 
relaxation time defined as r = /q°° dt. Figure 3(a) and (b) show r as a function of 

A^ for /3 = 3 with S" = 2,3,4 for the (a) cascade system and (b) loop system. For S > 3, 
the relaxation time r increases by several orders of magnitude with a decrease in A^ in both 
systems. On the other hand, r for = 2 does not exhibit any drastic change with the 
decrease in A^ in both systems. 

This prolongation of r for S > 2 becomes more prominent as (3 is increased. From several 
data, r is suggested to increase as a function of exp{f3V). Combining A^ and /3 dependencies. 
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we introduce a parameter h — N exp{—f3V). The discreteness effect is dominant when 
h — N exp{—/3V) is less than unity. Figure 3(c) and 3(d) show p — t/tn-^oo as a function of 
h for the (c) cascade system and (d) loop system for several values of f3 and 5" = 2, 3, 4. For 
S > 2, the deviation of p from the continuum limit (p = 1) becomes prominent when h is 
below unity in both systems. The increase in p appears to become steeper with an increase 
in S. On the other hand, p for 5" = 2 does not exhibit a drastic increase with a decrease in 
h. 



IV. ORIGIN OF SLOW RELAXATIONS AND CROSSOVER 
A. Relaxation processes for N ^ oo and N = 1 

Now, we analytically estimate the enhancement in relaxation time and explain its rep- 
resentation in the form h = Nexp{—pV). For this purpose, wc compare the estimate by 
Master equation analysis for small N and compare it with that from the continuum limit 
N ^ oo. 

In the continuum limit, the reaction dynamics are represented by the following rate 
equation: 

Xi^Xc[e'''^^{^-Xi) -Xi] (1) 

with Xi — pi/SN. Here, = 1 for i = 5* in the cascade system and Xc — xi for i — S 
in the loop system. In both systems, Xi — >■ x^'^ — s{i+e-pv) ^^^'^^ i — >■ oo. When the 
deviation from equilibrium 8xi small, its evolution for the loop systems obeys 

the following linearized equation 

Sxi = —Sxi. (2) 

For the cascade system, this equation is also valid for the elements i ^ S, whereas 5xs — 
—5xs- Then, the auto-correlation function of a small fluctuation of pj around is obtained 
as 



C{t) = exp(--^i) (3) 



for the loop system, and 



1 ^ _ 1 e-^^ 
C{t) = - exp(-i) + exp( —t) (4) 
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for the cascade system. Indeed, these agree quite well with the simulation results for a 
sufficiently large N (e.g., N — 1024 in Fig. 2.). Thus, the characteristic time of the 
relaxation is estimated as t^{S) ~ Se^^ for the loop system and r'^{S) ~ ^ + (S* — l)e^^ 
for the cascade system, which are consistent with the simulation results shown in Fig. 3. 

As the other extreme limit, consider the case with = 1. In this case, the relaxation 
dynamics are dominated by a completely different process induced by the absence of catalysts 
whose number can often go to zero. In such cases, states are trapped at some local energy 
minimum that appears due to the deficiency of catalysts. Then, the hopping processes 
among them play an important role in the relaxation dynamics, as shown below. In the 
following, we focus the cases with S* = 2 and S* = 3 to clarify that such an effect is induced 
by discreteness in the molecule number. Note that, as shown in the last section, the behavior 
for 5* > 3 is distinct from that for 5* = 2; in the former case, the relaxation time is enhanced 
by the decrease in A^, in contrast to the latter case. 

First, we study the loop system. When S — 2, the system reahzes 3 states from the 
initial conditions — (pi,P2) = (1,0), (0,1), and (1,1) — as shown in Fig. 4(a). Then, we 
estimate the time of the transition between (1, 0) and (0, 1). First, the transition rate from 
the state (1,0) to (1,1) is estimated as follows: for this transition, a pair of molecules from 
the product of the first species and the resource of the second species has to be chosen. This 
probability is given by §2^' while the reaction rate is given by e"^^ . Hence the rate is 
given by 2 • la^e"'^^ = e~^^ . Thus, the characteristic time of the correlation of each pi is 
given by ~ e^^ , which is consistent with the results shown in Fig. 2(d). 

On the other hand, for S = 3, the system realizes 7 states — (^1,^2,^3) = (1, 0, 0), (0, 1, 0), 
(0, 0, 1) ,(1, 1, 0), (1, 0, 1), (0, 1,1), and (1, 1, 1)— as shown in Fig. 4(b). The characteristic 
time of the correlation of each pi is given by the transition time among the three branches 
including lowest-energy states, (1, 0, 1) - (1, 0, 0), (1, 1, 0) - (0, 1, 0), and (0, 1, 1) - (0, 0, 1). 
Here, in order to hop from one branch to another, the system must go through the highest- 
energy state (1, 1, 1), due to the restriction by the catalytic relation. Now, we define the 
probability that the states in the branch (1,0,1) - (1,0,0) are realized as Qi,o,o- Then, 
the probability to realize the state (1,0,1) is given by ^^f-^fsv Q 1,0,0- Here, the transition 
rate from (1,0, 1) to (1, 1, 1) is given by Then, the probability current from the the 

branch (1,0,1) - (1,0,0) is estimated by - ^1^=^^1,0,0 ~ le-^^^Qi,o,o (e"'^^ « 1). 
Because of the symmetry among the catalytic reactions, the probability currents from the 
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FIG. 4: (a) Illustration of transition diagrams of (a) loop system with S* = 2, (b) loop system with 
S = s, (c) cascade system with 5 = 2, and (d) cascade system with 5 = 3, where arrows indicate 
possible transitions and the values next to them specify the transition ratios. 

other branches are obtained in the same way, to get the same form. Thus, the escape rate 
from each branch is estimated by ~ ^e~'^^^, and the characteristic time of the correlation 
of each pi is estimated as ~ 2e^^^. Because the relaxation time in the continuum limit is 
proportional to exp{l3V), the deviation p from it increases with exp{l3V), which is consistent 
with the results shown in Fig. 2(e). Thus, the enhancement of the relaxation time from the 
continuous case is explained. 

Essentially the same argument is also valid for the cascade systems. When S — 2, the 
system can realize transitions among 4 states — (0, 1) — (1, 1) — (1, 0) — (0, 0) — as shown in Fig. 
4(c). Here, (0, 1) is a metastable state and (0, 0) is the lowest-energy state. The relaxation 
is characterized by the escape rate from a metastable state, which is given by ~ e~^^ . Thus, 
the characteristic time of the correlation of each pi is given by ~ e^^' . 

On the other hand, for S = 3, the system realizes 8 states~(pi,p2,P3) = (0, 0, 0), (1, 0, 0), 
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(0, 1, 0), (0, 0, 1) ,(1, 1, 0), (1, 0, 1), (0, 1, 1), and (1, 1, 1)— as shown in Fig. 4(d). The slowest 
characteristic time of the relaxation is given by the transition time from the branch, (1, 1, 0) 
- (0,1,0) since the system must go through the highest-energy state (1,1,1), which is a 
limiting process for this case. Then, in a manner similar to the loop system with S = 3, the 
characteristic time is obtained as 2e^''^^\ This gives the characteristic time of the slowest 
motions of the system. This estimation fits well with the numerical result shown in Fig. 
2(b). 



B. N, j3 dependencies of C{oo) and relaxation time 

Next, we extend the argument of the last subsection to analyze the N and /3 dependencies 
of C(oo) and the relaxation time in greater detail. In particular, we explain why h = 
N exp{l3V) ~ 1 gives a critical value and how the amplification of relaxation time depends 
on /i for /i < 1. Because of the simplicity due to the symmetry in the catalytic relationship, 
we only study loop systems; however, the argument presented below can be extended to 
cascade systems. 

Figure 5(a) shows the transition diagram of the loop system with S = 2, where each 
circle indicates each state {pi,P2) and the arrows indicate possible transitions. Generally, 
for any values of S, the transition rate from a state {pi,P2, ■■■,Pi = n,pi+i, -.-yPs) to a state 
(pi,P2, ■■■,Pi — n + l,Pi+i, ■■■,Ps) per unit time is estimated as follows. For this transition, a 
pair of molecules from the resource of the iih species (Ri) and the product of the (i + l)th 
species (Pi+i) has to be chosen. This probabihty is given by sn-i ' reaction rate 

is given by e^^. Hence, the transition rate per unit time is given by l^/i^^i+i = ^'^gl^^^i'^^ e~^'^. 
Similarly, the transition rate in the opposite direction is given as = ^""g^^'i^ ■ If 

the molecule number is so large or /3 is so small that h = Ne~^^ >> 1, > W^ri+i-+n 

holds for small n and < holds for large n. Then, the dominant states of 

the system are located in an intermediate region in the phase space [0, A?"]. For example, the 
blue region in Fig. 5(a) indicates such dominant states for S — 2. 

Now, we define the probability that Pi — n as Q\^, and the joint probability to realize 
Pi = n and pi+i = m as Q\^. Here, Q\ = EZ=oQn,m and Qj,^^ = QnQm^- Then the time 
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evolution of (5n,m follows 

Q\m = J^^[{N-{n-l))e-^''(^^_,^^-r{n-rl)QU,,^-n^ (5) 
Then, we obtain 

where < pi >— Y^n=o''^Qn (< Pj+i Z]m=o"^Qm^)- Using this equation, we obtain the 
time evolution of < pj > as 

<Pi>- ^^[- <Pi> +iN- < Pi >)e-n (7) 

This implies that Xi =< pi > /SN obeys equation (1) for a sufficiently large value of N. 

On the other hand, if is so small or f3 is so large that h « 1, W^^^_^_i « 
holds for all i and n. Thus, Pi for all i tend to decrease to 0. Then, there exist 5"^ metastable 
states— (n, 0, 0, 0), (0, n, 0, 0), ... , (0, 0, 0, n, 0, .., 0), ... , and (0, 0, 0, n) {1 < n < 
N). Among them, the following 5" states, (1, 0, 0, ...0), (0, 1, 0, ...0), and (0, 0, 0, ...1), have 
the lowest energy. For example, in the cases with S — 2, the states (0,p) and (g, 0) (p, g 7^ 0) 
are metastable states and (1,0) and (0, 1) are the lowest-energy states. 

It should be noted that the lowest-energy states arc the dominant states for h << 1. 
The probability to realize these lowest-energy states tends to 1/S with an increase in /3. 
Thus, with the increase in /3, < Pi> approaches 1/S for small N, which indicates C(oo) = 
const. > for small N and large /3. 

Moreover, for /i << 1, the transitions among lowest-energy states contribute dominantly 
to the relaxation process. Then, we estimate the characteristic time of the fluctuations of 
the system for /i << 1 by considering the transition processes from one lowest-energy states 
such as (0, 0, 0,pj = 1,0...,0) to the other lowest-energy states such as (0, 0...,0,pj = 
0,0, pj/ = 1,0, ...,0). In the following, we consider only the cases with S — 2 and S — 3. 
We only focus on the dynamics of pj under the constraint that pj has only or 1, because 
h « 1. 

First, consider the case with S = 2. Figure 5(b) shows a detailed transition diagram 

around the region where pi (i = 1, 2) arc only or 1. The escape rate from (1, 0) and (0, 1) 
are given by ~ 2N-i ^~^'^ ■ Thus, the characteristic time of the correlation of each pi is given 

by 

r^(2) - ^e'^^ (8) 
12 
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FIG. 5: (a) Illustration of the transition diagrams for 5 = 2, and effective transition diagrams for 
(b) S = 2 and (c) 5 = 3, where bold arrows indicate the focused transitions in the text. 

which is consistent with the results shown in Fig. 6(a). 

Next, we study the case with 5* = 3. The transition diagram of the states (j5i,j925P3) 
is shown in Fig.5(c) when pi {i = 1,2,3) take only or 1. Similar to the = 1 case, 
the characteristic time of the transition among the three branches including lowest-energy 
states, (1, 0, 1) - (1, 0, 0), (1, 1, 0) - (0, 1, 0), and (0, 1, 1) - (0, 0, 1) through the state (1, 1, 1) is 
considered. In a manner similar to the = 1 case, the transition rate from each branch is es- 
timated by ~ i+jve-^v ~ (3N-i)(i+Ne-iiv) ' Thus, the relaxation time of the fluctuation 
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FIG. 6: Relaxation time r obtained from simulations (points) and its approximate analytical 
expression ^"^"(5) (curves) estimated in the text. Plotted as a function of N for the loop systems 
with (a) = 2 and (b) S = 3 with /3 = 2, 3, 4. The analytical expression agrees with the simulation 
data both for small N, and for large N, where r approaches a constant value expected from the 
rate equation. The crossover occurs at around h = Ne~^^ ~ 1 {N ~ exp(2), exp(3), and exp(4) 
for /? = 2, 3, and 4.). 

of pi is estimated as the decrease with as 

Considering the e^^ dependence of r^^oo, the above estimate is consistent with Fig. 6(b). 

For S larger than 3, the transition diagram becomes rather complicated. However, a 
similar analysis should be possible to estimate the prolongation in the relaxation time. 
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V. SUMMARY AND DISCUSSIONS 



In the present paper, the slowing down of the relaxation in reversible catalytic reaction 
networks induced by the smallness of molecule number is investigated as a general property 
of catalytic reaction networks. This prolongation of relaxation is a result of bottlenecks 
in reactions; these appear due to the deficiency of the catalyst required for a reaction. 
The number of molecules can be so small that the number of catalysts becomes zero. In 
this case, a pair of a substrate and the corresponding catalyst molecule species can hardly 
exist simultaneously. Such a constraint makes it difficult to realize a specific configuration 
necessary for the relaxation. The probability for realization is given by exp{—PEhottie), with 
Ebottie as the corresponding energy barrier to realize such rare conditions, or the sum of 
such energy barriers. This bottleneck energy is generally different from the energy gap in 
the continuum limit that is obtained from the rate equation (ordinary differential equation). 
Thus, the relaxation time at a small molecule number deviates from the continuum case by 
the factor exp{l36E) with an appropriate effective energy difference, 6E. 

By considering the models of simple catalytic reaction networks consisting of resource 
chemicals of S species and the corresponding products, we have demonstrated this devia- 
tion of relaxation time from both direct simulations and analysis by using Master equation. 
From the numerical and analytic estimates, Ef,ouieneck = '^V and 6E = V for S = 3, where 
V is the energy gap between the resource and the product chemicals. For S > 2, in general, 
the prolongation of the relaxation time becomes prominent when h = N exp{—l3V) is less 
than unity, and its amplification ratio from the continuum limit is represented as a function 
of 5* and h. Note that the cascade system in the = 1 case is equivalent to the "Asymmet- 
rically Constrained Ising Chain" (ACIC), Hierarchically constrained Ising model, or East 
model, which are studied as simple abstract models for glassy states 23l-l25l|. Following the 
interpretation therein, the increase in relaxation time at < 1 as a result of the decrease 
in A^ or temperature may be regarded as a type of glass transition. According to the recent 
studies on ACIC, the correlation time of the motion of pi (not the relaxation time of the 
;otal _system) is estimated as ri ~ (1 + 2e^^)'^ where the integer k obeys 2*^^^^) < S < 2^ 
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25j. In cases with = 2, 3, 4, this fact is consistent with our estimate of the relaxation 



time of the cascade system with A^ = 1. The estimation of 6E as a function of S and h for 
general cases both for cascade and loop systems is an important issue that should be studied 
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in the future. 

In addition to the slow down in relaxation, the equilibrium distribution deviates in a 
network called a loop system, where all the reactions are catalyzed by one of the products. 
The constraint that the numbers of a certain pair of chemical species cannot simultaneously 
be zero leads to the deviation of the average distribution of molecule numbers from the 
continuum limit. Again, this deviation becomes prominent when h is less than unity. 

Although we have adopted simple network motives to analyze the relaxation, the pro- 
longation of relaxation time is quite general in catalytic reaction networks. Catalytic bot- 
tlenecks often appear as the number of molecules is decreased in a large variety of reaction 



networks in which catalysts are synthesized within |21l. |22|. The present study can provide a 
basis for the general case with complex networks, as the motives here are sufficiently small 
and can exist within such complex networks. 

Biochemical reactions generally progress in the presence of catalysts that are themselves 
synthesized as products of such reactions. These reactions form a network of a variety of 
chemical species. Here, the molecule number of each species is generally not very large. 
Hence, the slow relaxation process and deviation from equilibrium discussed in this study 
may underlie intracellular reaction processes. Moreover, the present network motives are 
so simple that they are suggested to exist in biochemical networks. We also note that the 
resource and product in our model can be interpreted as non-excited and excited states 
of enzymatic molecules. Indeed, many molecules are known to exhibit catalytic activity 
only when they are in an excited state, which can help other chemicals to switch to an 
excited state. In fact, such networks with mutual excitation are known in signal-transduction 
networks 26|-|28| , where the present slow relaxation mechanism may be relevant to sustain 



the excitability of a specific enzyme type over a long time span. It is important to pursue the 
relevance of the present mechanism in cell-biological problems by considering more realistic 
models in the future. 

We also note that not only the discreteness in the molecule number but also the negative 
correlation between a substrate and the corresponding catalyst within a reaction network 



22j. The 



or in a spatial concentration pattern suppresses the relaxation process 
present mechanism due to discreteness may work synergetically with the earlier mechanism 
to further suppress the relaxation to equilibrium. The construction of reaction networks to 
achieve slower relaxation together with the network analysis will be an important issue in 
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the future. 
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